
% skymaps.m
% compute and plot antenna response patterns for different detectors,
% using the coherent network code developed at Caltech in 2005
% in lscsoft/src/matapps/src/searches/burst/coherent-network

% Tell matlab where the coherent-network code is 
%  (or you can do this with the MATLABPATH environment variable)
addpath /archive/home/ajw/lscsoft/src/matapps/src/searches/burst/coherent-network -begin
% make sure code in our local directory is picked up first
addpath . -begin

% make a source sky grid with good resolution:
angularResolution = 0.05;
skyPositions = sinusoidalMap(angularResolution);

% extract the spherical angles associated with the sky grid:
thetas = skyPositions(:,1);
phis = skyPositions(:,2);
% we will average over plus and cross polarizations,
%  so the result doesn't depend on the polarization of the source.
psi = 0;

% define the map structure:
projectionType = 'mollweid';
mapstruct = defaultm(projectionType);
mapstruct.angleunits = 'degrees';
mapstruct.origin = [0 0 0];
mapstruct.falseeasting = 0;
mapstruct.falsenorthing = 0;
mapstruct.scalefactor = 1;

% load the worldmap coastlines:
load coast;
[xc, yc] = mfwdtran(mapstruct, lat, long);

% define a set of detectors recognized by LoadDetectorData
detectors{1} = 'LHO';
detectors{2} = 'LLO';
detectors{3} = 'Virgo';
%detectors{4} = 'INDIGO';
Ndetectors = length(detectors);

% you can define your own detectors.
% For example, let's put INDIGO in Kathmandu:
%     Det.h = 1355.;           % elevation in meters
%     Det.phi = 27.699346;     % latitude in degrees
%     Det.lambda = 85.316831;  % longitude in degrees
%     thet = 0.5;              % angle of Y arm wrt west, in radians
% 
% % construct the detector definition from the above numbers:
%     deg2rad = pi/180.;
%     Rearth = 6.378e6;
%     RR = Rearth+Det.h;
%     %------ Location of vertex, in Cartesian Earth-based coordinates (m)
%     DetX = RR*cos(deg2rad*Det.phi)*cos(deg2rad*Det.lambda);
%     DetY = RR*cos(deg2rad*Det.phi)*sin(deg2rad*Det.lambda);
%     DetZ = RR*sin(deg2rad*Det.phi);
%     Det.V = [ DetX ; DetY ; DetZ ];
%     %----- zenith vector.  
%     Det.Z = Det.V/norm(Det.V);
%     %----- Arm pointing directions
%     NorthPole = [0; 0; 1];
%     Det.Y = cross(Det.Z,NorthPole);
%     Det.X = cross(Det.Y,Det.Z);
%     Det.X = Det.X/norm(Det.X);
%     Det.Y = Det.Y/norm(Det.Y);
%     %----- Rotate by some angle from due west:
%     DetX = cos(thet)*Det.X+sin(thet)*Det.Y;
%     DetY =-sin(thet)*Det.X+cos(thet)*Det.Y;
%     Det.X = DetX/norm(DetX);
%     Det.Y = DetY/norm(DetY);
%     %----- Detector matrix:
%     Det.d = 0.5*(Det.X*Det.X' - Det.Y*Det.Y') ;
%     Det.name = 'INDIGO';
%     Det.n = [];
%     Det.psi = [];
%     %----- Save this detector definition
%     INDIGO = Det;


% loop over detector
for I=1:Ndetectors,
   
%    if (strcmp(detectors{I},'INDIGO'))
%      Det = INDIGO;
%    else
     Det = LoadDetectorData(detectors{I});
%    end
    

% compute antenna response
   [Fp, Fc] = ComputeAntennaResponse(phis,thetas,psi,Det.d);
% polarization-averaged response:
   FF = sqrt(Fp.^2+Fc.^2);
   
   if (I==1)
       R = 0;
   else
       R = R + FF.^2;
   end
   
   

% and plot it with xproject:
   figure(I)
   handle1 = xproject(skyPositions, FF, projectionType);
   title(detectors{I});
   colorbar
   

% add world map and a dot for the location of the detector:
  [x, y] = mfwdtran(mapstruct, Det.phi, Det.lambda);
  hold on
  plot(x,y,'k*');
  plot(xc,yc);
  hold off
   end

   figure(12)
   handle1 = xproject(skyPositions, sqrt(R), projectionType);
   title('Network');
   colorbar
   

% add world map and a dot for the location of the detector:
  [x, y] = mfwdtran(mapstruct, Det.phi, Det.lambda);
  hold on
  plot(x,y,'k*');
  plot(xc,yc);
  hold off
